CpG Methylation Analysis using M3D

M3D is an R package that distinguishes DMR (differential methylation regions) by kernel methods. (SVM?) Promoter regions of housekeeping genes and its CpG islands are known to be less methylated. This makes a difference in gene expression regulation. Therefore, promoter region’s degree of methylation is an indicator of gene expression.

Also, RRBS (Reduced Representation Bisulfite Sequencing) is an efficient experiment for measuring methylation of dense CpG regions’ cytosines.

  1. Digestion: MspI is used to digest genomic DNA, unbiased towards CpG methylation. It targets 5’CCGG3’ sequences and cleaves phosphodiester bonds upstream of CpG dinucleotide. Each fragment has a CpG at the end then.

  2. End repair and A-tailing: MspI cleaves double stranded DNA and results in sticky ends. End-repair is required. And then add A-tail to both plus and minus strands at the 3’ terminal. It’s necessary for adapter ligation. (dATPs are added in excess to increase efficiency.)

  3. Sequence adapters: Methylated sequence adapters (oligonucleotides) -> All cytosines are replaced with 5’methyl-cytosines.

  4. Fragment purification: Gel electrophoresis and select for 40-220 bp length fragments that are representiative of the majority of promoter sequences and CpG islands.

  5. Bisulfite conversion: Application of bisulfite results in the deamination of unmethylated cytosines into uracils.

  6. PCR amplification

  7. PCR purification

  8. Sequencing

Statistical analysis includes β binomial distribution hypothesis testing for CpG loci, and individual CpGs are connected in order to produce an estimated values of differentially methylated regions (DMR).

Install Bioconductor Manager, M3D, BiSeq

my_packages <- c("M3D", "BiSeq")
lapply(my_packages, library, character.only=TRUE)
[[1]]
 [1] "BiSeq"                "Formula"              "SummarizedExperiment" "DelayedArray"        
 [5] "BiocParallel"         "matrixStats"          "Biobase"              "GenomicRanges"       
 [9] "GenomeInfoDb"         "IRanges"              "S4Vectors"            "BiocGenerics"        
[13] "parallel"             "stats4"               "M3D"                  "stats"               
[17] "graphics"             "grDevices"            "utils"                "datasets"            
[21] "methods"              "base"                

[[2]]
 [1] "BiSeq"                "Formula"              "SummarizedExperiment" "DelayedArray"        
 [5] "BiocParallel"         "matrixStats"          "Biobase"              "GenomicRanges"       
 [9] "GenomeInfoDb"         "IRanges"              "S4Vectors"            "BiocGenerics"        
[13] "parallel"             "stats4"               "M3D"                  "stats"               
[17] "graphics"             "grDevices"            "utils"                "datasets"            
[21] "methods"              "base"                
group <- factor(c('H1-hESC', 'H1-hESC', 'K562', 'K562'))
samples <- c('H1-hESC1', 'H1-hESC2', 'K562-1', 'K562-2')
colData <- data.frame(group, row.names = samples)
colData
ABCDEFGHIJ0123456789
 
 
group
<fctr>
H1-hESC1H1-hESC
H1-hESC2H1-hESC
K562-1K562
K562-2K562
4 rows
levels(colData$group)
[1] "H1-hESC" "K562"   
class(colData)
[1] "data.frame"

Note that data.frame and DataFrame are different. data.frame is a built-in data structure, and DataFrame is a function in S4Vectors and utilized by Bioconductor.

colData <- DataFrame(group, row.names = samples)
class(colData)
[1] "DFrame"
attr(,"package")
[1] "S4Vectors"

DataFrame is a data structure included in Bioconductor. Refer to this link: Bioconductor DataFrame explanation

The file that will be used in this session are ENCODE RRBS data of H1-hESC and K562 cell lines. H1-hESC cell’s two samples are two be compared to two samples of K562 leukemia cell. Download the necessary bed file from UCSC (http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/)

A BED file (.bed) is a tab-delimited text file that defines a feature track. It can have any file extension, but .bed is recommended. The BED file format is described on the UCSC Genome Bioinformatics web site: http://genome.ucsc.edu/FAQ/FAQformat. Tracks in the UCSC Genome Browser (http://genome.ucsc.edu/) can be downloaded to BED files and loaded into IGV.

path = "/home/seungho/GDA/PublicENCODE" # path to the directory
setwd(path)
The working directory was changed to /home/seungho/GDA/PublicENCODE inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.
files <- c('wgEncodeHaibMethylRrbsH1hescHaibSitesRep1.bed.gz', 'wgEncodeHaibMethylRrbsH1hescHaibSitesRep2.bed.gz', 'wgEncodeHaibMethylRrbsK562HaibSitesRep1.bed.gz', 'wgEncodeHaibMethylRrbsK562HaibSitesRep2.bed.gz')
rrbs <- readBedFiles(files, colData, bed_type = 'Encode')
Reading in files
Processing sample H1-hESC1 ... 
Read 1288079 records
Begin sorting data
Data sort complete
Processing sample H1-hESC2 ... 
Read 1302577 records
Begin sorting data
Data sort complete
Processing sample K562-1 ... 
Read 1275818 records
Begin sorting data
Data sort complete
Processing sample K562-2 ... 
Read 1114927 records
Begin sorting data
Data sort complete
Building BSraw object.

BED format file is, again, a file format that you are able to download through UCSC’s ENCODE data portal. It is tab-delimited and includes next attributes:

the chromosome name, start coordinate, end coordinate, some identifier, read coverage (#T), strand, start and end coordinates again(we discard this duplicated information), some color value, read coverage(#T) and the methylation percentage.

BED format file must have a header row, and the coordinate reference for the locus is 0.

readBedFile is a function that resides in the namespace of M3D. The following is its code.

function (files, colData, bed_type = "Encode", eData = NaN) 
{
    if (nrow(colData) != length(files)) {
        stop("Row number of colData must equal length of files.")
    }
    if (!test_colData_structure(colData)) {
        stop("Files must be specified group by group. Please re-make colData to\n             show all of one group, then all of the next, etc.")
    }
    message("Reading in files\n")
    methData <- readBedRaw(files, colData, bed_type = bed_type)
    cat("Building BSraw object.\n")
    fData <- methData[[1]]
    if (length(methData) > 1) {
        for (i in seq(along = methData)[-1]) {
            fData <- unique(c(fData, methData[[i]]))
        }
    }
    elementMetadata(fData) <- NULL
    names(fData) <- as.character(1:length(fData))
    tReads <- matrix(integer(length = length(fData) * length(methData)), 
        nrow = length(fData))
    mReads <- matrix(integer(length = length(fData) * length(methData)), 
        nrow = length(fData))
    for (i in seq(along = methData)) {
        mtch <- findOverlaps(fData, methData[[i]])
        mtch.m <- as.matrix(mtch)
        ind1 <- mtch.m[, 1]
        ind2 <- mtch.m[, 2]
        tReads[ind1, i] <- elementMetadata(methData[[i]])$reads[ind2]
        mReads[ind1, i] <- elementMetadata(methData[[i]])$methylated[ind2]
    }
    colnames(tReads) <- rownames(colData)
    colnames(mReads) <- rownames(colData)
    rownames(tReads) <- names(fData)
    rownames(mReads) <- names(fData)
    if (is.na(eData)) {
        eData <- SimpleList(exp = "Demo")
    }
    rrbs = BSraw(exptData = eData, colData = colData, rowRanges = fData, 
        totalReads = tReads, methReads = mReads)
    return(rrbs)
}
<bytecode: 0x55df317be410>
<environment: namespace:M3D>

Next, use BiSeq package to generate GRanges file for the region. BiSeq package is built for processing bisulfite-conversion sequenced data.

rrbs.CpGs <- clusterSites(object = rrbs, groups = unique(colData(rrbs)$group), perc.samples = 3/4, min.sites = 20, max.dist = 100)
CpGs <- clusterSitesToGR(rrbs.CpGs)

In order to understand the next step, you gotta know what coverage (or read depth) is. Coverage in DNA sequencing is the number of unique reads that include a given nucleotide in the reconstructed sequence. Deep sequencing aims at higher coverage of each unique reads.

Now, we cut off the regions with coverage (read depth) less than 100. totalReads(rrbs.CpGs) is a matrix that contains position and total coverage of each position.

In the following code, overlaps is

> overlaps
Hits object with 576106 hits and 0 metadata columns:
           queryHits subjectHits
           <integer>   <integer>
       [1]         1           1
       [2]         1           2
       [3]         1           3
       [4]         1           4
       [5]         1           5
       ...       ...         ...
  [576102]     14453      576102
  [576103]     14453      576103
  [576104]     14453      576104
  [576105]     14453      576105
  [576106]     14453      576106
  -------

and

subjectHits(overlaps[ queryHits(overlaps)== 1])
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

subjectHits brings indices of where that number of overlap is.

inds <- vector()
overlaps <- findOverlaps(CpGs, rowRanges(rrbs.CpGs))
for (i in 1:length(CpGs)){
  covs <- colSums(totalReads(rrbs.CpGs)[subjectHits(overlaps[ queryHits(overlaps) == i]),])
  if (!any(covs <= 100)){
    inds <- c(inds, i)
  }
}
CpGs <- CpGs[inds]
rm(inds)

For next analysis, use only 1000 regions. Update overlaps object and explain correct overlaps.

rrbs <- rrbs.CpGs
CpGs <- CpGs[1:1000]
overlaps <- findOverlaps(CpGs, rowRanges(rrbs))
rrbs <- rrbs[subjectHits(overlaps)]
overlaps <- findOverlaps(CpGs, rowRanges(rrbs))

Use the two data in ENCODE consortium H1-hESC and K562 as the two RRBS files.

Use BiSeq package’s functions clusterSites and clusterSitesToGR for clustering.

Delete islands of all samples that have less than 100? and include only the first thousand.

Analysis is stored in GRanges and each item shows a region’s start and end.

rrbs
class: BSraw 
dim: 39907 4 
metadata(0):
assays(2): totalReads methReads
rownames(39907): 643464 643465 ... 687062 687063
rowData names(1): cluster.id
colnames(4): H1-hESC1 H1-hESC2 K562-1 K562-2
colData names(1): group
CpGs
GRanges object with 1000 ranges and 1 metadata column:
         seqnames              ranges strand | cluster.id
            <Rle>           <IRanges>  <Rle> |   <factor>
     [1]     chr1       840131-840357      * |     chr1_1
     [2]     chr1       845246-845561      * |     chr1_2
     [3]     chr1       858978-859375      * |     chr1_3
     [4]     chr1       859529-859727      * |     chr1_4
     [5]     chr1       875730-875953      * |     chr1_5
     ...      ...                 ...    ... .        ...
   [996]     chr1 168105869-168106338      * |  chr1_1009
   [997]     chr1 168148116-168148585      * |  chr1_1010
   [998]     chr1 168194959-168195245      * |  chr1_1011
   [999]     chr1 169075431-169075641      * |  chr1_1012
  [1000]     chr1 170633423-170633656      * |  chr1_1013
  -------
  seqinfo: 25 sequences from an unspecified genome; no seqlengths

Maximum Mean Discrepancy (MMD) and Coverage

M3D analysis uses MMD values for the entire methylation data and the applied regions. Both can be obtained through M3D_Wrapper function. Sample pairs’ entire methylated MMD matrix is the return value of M3D_Wrapper. The second result is equivalent coverage. MMDlist$Full and MMDlist$Coverage variables are those. Substitute the difference between the two values into M3Dstat, and this shows structure of each items.

head(MMDlist$Full)
     H1-hESC1  vs  H1-hESC2 H1-hESC1  vs  K562-1
[1,]             0.13328796            0.3711302
[2,]             0.15843597            0.1991325
[3,]             0.11149443            0.1376077
[4,]             0.05180868            0.1450703
[5,]             0.17830486            0.2438087
[6,]             0.02382480            0.0972950
     H1-hESC1  vs  K562-2 H1-hESC2  vs  K562-1 H1-hESC2  vs  K562-2
[1,]            0.3882027           0.23581746           0.22551184
[2,]            0.1901730           0.14194300           0.14213741
[3,]            0.1476828           0.03963481           0.05186073
[4,]            0.1196117           0.11307940           0.07325263
[5,]            0.1983181           0.07915218           0.10918123
[6,]            0.1378252           0.09148931           0.13047384
     K562-1  vs  K562-2
[1,]         0.06755053
[2,]         0.04260434
[3,]         0.01802223
[4,]         0.05798640
[5,]         0.09577036
[6,]         0.05080778
head(MMDlist$Coverage)
     H1-hESC1  vs  H1-hESC2 H1-hESC1  vs  K562-1 H1-hESC1  vs  K562-2 H1-hESC2  vs  K562-1 H1-hESC2  vs  K562-2 K562-1  vs  K562-2
[1,]             0.13245172            0.3979882            0.3998015           0.26189520           0.23546355         0.06338745
[2,]             0.21273851            0.1709892            0.1664227           0.05580718           0.07902125         0.04269711
[3,]             0.11562054            0.1408935            0.1493031           0.04014058           0.04976709         0.01636765
[4,]             0.05353733            0.1461974            0.1205286           0.11230110           0.07220178         0.05731698
[5,]             0.17977249            0.2456640            0.2034446           0.07999092           0.11754403         0.09997088
[6,]             0.02360063            0.0988395            0.1365794           0.09302259           0.12916722         0.04774588
M3Dstat <- MMDlist$Full - MMDlist$Coverage

Below is a scatterplot for comparison among replicates of Full MMD and Coverage MMD.

repCols <- c(1, 6)
plot(as.vector(MMDlist$Full[, repCols]), as.vector(MMDlist$Coverage[, repCols]), xlab = 'Full MMD', ylab = 'Coverage MMD', main = 'Test Statistic: Replicate Comparison')
abline(a = 0, b = 1, col='red', lwd=2)

groupCols <- 2:5
plot(as.vector(MMDlist$Full[, groupCols]), as.vector(MMDlist$Coverage[, groupCols]), xlab = 'Full MMD', ylab = 'Coverage MMD', main = 'Test Statistics: Inter-Group Comparison')
abline(a=0, b=1, col='red', lwd=2)

The following code prints M3D analysis result as histograms.

repCols <- c(1:6)
groupCols <- 2:5
M3Dstat <- MMDlist$Full - MMDlist$Coverage
breaks <- seq(-0.2, 1.2, 0.1)
hReps <- hist(M3Dstat[, repCols], breaks = breaks, plot = FALSE)
hGroups <- hist(rowMeans(M3Dstat[, groupCols]), breaks = breaks, plot = FALSE)
col1 <- "#0000FF40"
col2 <- "#FF000040"
plot(hReps, col = col1, freq = FALSE, ylab = "Density", xlab = "MMD statistic", main = "M3D Stat Histogram")
plot(hGroups, col = col2, freq = FALSE, add = TRUE)
legend(x = 0.5, y = 3, legend = c("Replicates", "Groups"), fill = c(col1, col2))

Next is histone modification data search in the EpiFactors database.

LS0tCnRpdGxlOiAiTTNEX3ByYWN0aWNlIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCiMgQ3BHIE1ldGh5bGF0aW9uIEFuYWx5c2lzIHVzaW5nIE0zRApNM0QgaXMgYW4gUiBwYWNrYWdlIHRoYXQgZGlzdGluZ3Vpc2hlcyBETVIgKGRpZmZlcmVudGlhbCBtZXRoeWxhdGlvbiByZWdpb25zKSBieSBrZXJuZWwgbWV0aG9kcy4gKFNWTT8pClByb21vdGVyIHJlZ2lvbnMgb2YgaG91c2VrZWVwaW5nIGdlbmVzIGFuZCBpdHMgQ3BHIGlzbGFuZHMgYXJlIGtub3duIHRvIGJlIGxlc3MgbWV0aHlsYXRlZC4gVGhpcyBtYWtlcyBhIGRpZmZlcmVuY2UgaW4gZ2VuZSBleHByZXNzaW9uIHJlZ3VsYXRpb24uIFRoZXJlZm9yZSwgcHJvbW90ZXIgcmVnaW9uJ3MgZGVncmVlIG9mIG1ldGh5bGF0aW9uIGlzIGFuIGluZGljYXRvciBvZiBnZW5lIGV4cHJlc3Npb24uCgpBbHNvLCBSUkJTIChSZWR1Y2VkIFJlcHJlc2VudGF0aW9uIEJpc3VsZml0ZSBTZXF1ZW5jaW5nKSBpcyBhbiBlZmZpY2llbnQgZXhwZXJpbWVudCBmb3IgbWVhc3VyaW5nIG1ldGh5bGF0aW9uIG9mIGRlbnNlIENwRyByZWdpb25zJyBjeXRvc2luZXMuCgoxLiBEaWdlc3Rpb246IE1zcEkgaXMgdXNlZCB0byBkaWdlc3QgZ2Vub21pYyBETkEsIHVuYmlhc2VkIHRvd2FyZHMgQ3BHIG1ldGh5bGF0aW9uLiBJdCB0YXJnZXRzIDUnQ0NHRzMnIHNlcXVlbmNlcyBhbmQgY2xlYXZlcyBwaG9zcGhvZGllc3RlciBib25kcyB1cHN0cmVhbSBvZiBDcEcgZGludWNsZW90aWRlLiBFYWNoIGZyYWdtZW50IGhhcyBhIENwRyBhdCB0aGUgZW5kIHRoZW4uCgoyLiBFbmQgcmVwYWlyIGFuZCBBLXRhaWxpbmc6IE1zcEkgY2xlYXZlcyBkb3VibGUgc3RyYW5kZWQgRE5BIGFuZCByZXN1bHRzIGluIHN0aWNreSBlbmRzLiBFbmQtcmVwYWlyIGlzIHJlcXVpcmVkLiBBbmQgdGhlbiBhZGQgQS10YWlsIHRvIGJvdGggcGx1cyBhbmQgbWludXMgc3RyYW5kcyBhdCB0aGUgMycgdGVybWluYWwuIEl0J3MgbmVjZXNzYXJ5IGZvciBhZGFwdGVyIGxpZ2F0aW9uLiAoZEFUUHMgYXJlIGFkZGVkIGluIGV4Y2VzcyB0byBpbmNyZWFzZSBlZmZpY2llbmN5LikKCjMuIFNlcXVlbmNlIGFkYXB0ZXJzOiBNZXRoeWxhdGVkIHNlcXVlbmNlIGFkYXB0ZXJzIChvbGlnb251Y2xlb3RpZGVzKSAtPiBBbGwgY3l0b3NpbmVzIGFyZSByZXBsYWNlZCB3aXRoIDUnbWV0aHlsLWN5dG9zaW5lcy4KCjQuIEZyYWdtZW50IHB1cmlmaWNhdGlvbjogR2VsIGVsZWN0cm9waG9yZXNpcyBhbmQgc2VsZWN0IGZvciA0MC0yMjAgYnAgbGVuZ3RoIGZyYWdtZW50cyB0aGF0IGFyZSByZXByZXNlbnRpYXRpdmUgb2YgdGhlIG1ham9yaXR5IG9mIHByb21vdGVyIHNlcXVlbmNlcyBhbmQgQ3BHIGlzbGFuZHMuCgo1LiBCaXN1bGZpdGUgY29udmVyc2lvbjogQXBwbGljYXRpb24gb2YgYmlzdWxmaXRlIHJlc3VsdHMgaW4gdGhlIGRlYW1pbmF0aW9uIG9mIHVubWV0aHlsYXRlZCBjeXRvc2luZXMgaW50byB1cmFjaWxzLgoKNi4gUENSIGFtcGxpZmljYXRpb24KCjcuIFBDUiBwdXJpZmljYXRpb24KCjguIFNlcXVlbmNpbmcKCiFbXShodHRwczovL3VwbG9hZC53aWtpbWVkaWEub3JnL3dpa2lwZWRpYS9jb21tb25zL2YvZmYvUmVkdWNlZF9SZXByZXNlbnRhdGlvbl9CaXN1bGZpdGVfU2VxdWVuY2luZ19Qcm90b2NvbC5qcGcpCiAKU3RhdGlzdGljYWwgYW5hbHlzaXMgaW5jbHVkZXMgJFxiZXRhJCBiaW5vbWlhbCBkaXN0cmlidXRpb24gaHlwb3RoZXNpcyB0ZXN0aW5nIGZvciBDcEcgbG9jaSwgYW5kIGluZGl2aWR1YWwgQ3BHcyBhcmUgY29ubmVjdGVkIGluIG9yZGVyIHRvIHByb2R1Y2UgYW4gZXN0aW1hdGVkIHZhbHVlcyBvZiBkaWZmZXJlbnRpYWxseSBtZXRoeWxhdGVkIHJlZ2lvbnMgKERNUikuCgojIyBJbnN0YWxsIEJpb2NvbmR1Y3RvciBNYW5hZ2VyLCBNM0QsIEJpU2VxCmBgYHtyLCBlY2hvPUZBTFNFfQppbnN0YWxsLnBhY2thZ2VzKCJCaW9jTWFuYWdlciIpCkJpb2NNYW5hZ2VyOjppbnN0YWxsKCJNM0QiKQpCaW9jTWFuYWdlcjo6aW5zdGFsbCgiQmlTZXEiKQpgYGAKCgpgYGB7cn0KbXlfcGFja2FnZXMgPC0gYygiTTNEIiwgIkJpU2VxIikKbGFwcGx5KG15X3BhY2thZ2VzLCBsaWJyYXJ5LCBjaGFyYWN0ZXIub25seT1UUlVFKQpgYGAKYGBge3J9Cmdyb3VwIDwtIGZhY3RvcihjKCdIMS1oRVNDJywgJ0gxLWhFU0MnLCAnSzU2MicsICdLNTYyJykpCnNhbXBsZXMgPC0gYygnSDEtaEVTQzEnLCAnSDEtaEVTQzInLCAnSzU2Mi0xJywgJ0s1NjItMicpCmNvbERhdGEgPC0gZGF0YS5mcmFtZShncm91cCwgcm93Lm5hbWVzID0gc2FtcGxlcykKY29sRGF0YQpgYGAKCmBgYHtyfQpsZXZlbHMoY29sRGF0YSRncm91cCkKYGBgCmBgYHtyfQpjbGFzcyhjb2xEYXRhKQpgYGAKTm90ZSB0aGF0IGBkYXRhLmZyYW1lYCBhbmQgYERhdGFGcmFtZWAgYXJlIGRpZmZlcmVudC4gYGRhdGEuZnJhbWVgIGlzIGEgYnVpbHQtaW4gZGF0YSBzdHJ1Y3R1cmUsIGFuZCBgRGF0YUZyYW1lYCBpcyBhIGZ1bmN0aW9uIGluIFM0VmVjdG9ycyBhbmQgdXRpbGl6ZWQgYnkgQmlvY29uZHVjdG9yLgoKYGBge3J9CmNvbERhdGEgPC0gRGF0YUZyYW1lKGdyb3VwLCByb3cubmFtZXMgPSBzYW1wbGVzKQpjbGFzcyhjb2xEYXRhKQpgYGAKRGF0YUZyYW1lIGlzIGEgZGF0YSBzdHJ1Y3R1cmUgaW5jbHVkZWQgaW4gQmlvY29uZHVjdG9yLiBSZWZlciB0byB0aGlzIGxpbms6CltCaW9jb25kdWN0b3IgRGF0YUZyYW1lIGV4cGxhbmF0aW9uXShodHRwczovL3d3dy5iaW9jb25kdWN0b3Iub3JnL2hlbHAvY291cnNlLW1hdGVyaWFscy8yMDE5L0Jpb2NEZXZlbEZvcnVtLzAyLURhdGFGcmFtZS5wZGYpCgpUaGUgZmlsZSB0aGF0IHdpbGwgYmUgdXNlZCBpbiB0aGlzIHNlc3Npb24gYXJlIEVOQ09ERSBSUkJTIGRhdGEgb2YgSDEtaEVTQyBhbmQgSzU2MiBjZWxsIGxpbmVzLiBIMS1oRVNDIGNlbGwncyB0d28gc2FtcGxlcyBhcmUgdHdvIGJlIGNvbXBhcmVkIHRvIHR3byBzYW1wbGVzIG9mIEs1NjIgbGV1a2VtaWEgY2VsbC4gRG93bmxvYWQgdGhlIG5lY2Vzc2FyeSBiZWQgZmlsZSBmcm9tIFVDU0MgKGh0dHA6Ly9oZ2Rvd25sb2FkLmNzZS51Y3NjLmVkdS9nb2xkZW5QYXRoL2hnMTkvZW5jb2RlRENDL3dnRW5jb2RlSGFpYk1ldGh5bFJyYnMvKQoKQSBCRUQgZmlsZSAoLmJlZCkgaXMgYSB0YWItZGVsaW1pdGVkIHRleHQgZmlsZSB0aGF0IGRlZmluZXMgYSBmZWF0dXJlIHRyYWNrLiBJdCBjYW4gaGF2ZSBhbnkgZmlsZSBleHRlbnNpb24sIGJ1dCAuYmVkIGlzIHJlY29tbWVuZGVkLiBUaGUgQkVEIGZpbGUgZm9ybWF0IGlzIGRlc2NyaWJlZCBvbiB0aGUgVUNTQyBHZW5vbWUgQmlvaW5mb3JtYXRpY3Mgd2ViIHNpdGU6IGh0dHA6Ly9nZW5vbWUudWNzYy5lZHUvRkFRL0ZBUWZvcm1hdC4gVHJhY2tzIGluIHRoZSBVQ1NDIEdlbm9tZSBCcm93c2VyIChodHRwOi8vZ2Vub21lLnVjc2MuZWR1LykgY2FuIGJlIGRvd25sb2FkZWQgdG8gQkVEIGZpbGVzIGFuZCBsb2FkZWQgaW50byBJR1YuCgpgYGB7cn0KcGF0aCA9ICIvaG9tZS9zZXVuZ2hvL0dEQS9QdWJsaWNFTkNPREUiICMgcGF0aCB0byB0aGUgZGlyZWN0b3J5CnNldHdkKHBhdGgpCmZpbGVzIDwtIGMoJ3dnRW5jb2RlSGFpYk1ldGh5bFJyYnNIMWhlc2NIYWliU2l0ZXNSZXAxLmJlZC5neicsICd3Z0VuY29kZUhhaWJNZXRoeWxScmJzSDFoZXNjSGFpYlNpdGVzUmVwMi5iZWQuZ3onLCAnd2dFbmNvZGVIYWliTWV0aHlsUnJic0s1NjJIYWliU2l0ZXNSZXAxLmJlZC5neicsICd3Z0VuY29kZUhhaWJNZXRoeWxScmJzSzU2MkhhaWJTaXRlc1JlcDIuYmVkLmd6JykKcnJicyA8LSByZWFkQmVkRmlsZXMoZmlsZXMsIGNvbERhdGEsIGJlZF90eXBlID0gJ0VuY29kZScpCmBgYApCRUQgZm9ybWF0IGZpbGUgaXMsIGFnYWluLCBhIGZpbGUgZm9ybWF0IHRoYXQgeW91IGFyZSBhYmxlIHRvIGRvd25sb2FkIHRocm91Z2ggVUNTQydzIEVOQ09ERSBkYXRhIHBvcnRhbC4gSXQgaXMgdGFiLWRlbGltaXRlZCBhbmQgaW5jbHVkZXMgbmV4dCBhdHRyaWJ1dGVzOgoKdGhlIGNocm9tb3NvbWUgbmFtZSwgc3RhcnQgY29vcmRpbmF0ZSwgZW5kIGNvb3JkaW5hdGUsIHNvbWUgaWRlbnRpZmllciwgcmVhZCBjb3ZlcmFnZSAoI1QpLCBzdHJhbmQsIHN0YXJ0IGFuZCBlbmQgY29vcmRpbmF0ZXMgYWdhaW4od2UgZGlzY2FyZCB0aGlzIGR1cGxpY2F0ZWQgaW5mb3JtYXRpb24pLCBzb21lIGNvbG9yIHZhbHVlLCByZWFkIGNvdmVyYWdlKCNUKSBhbmQgdGhlIG1ldGh5bGF0aW9uIHBlcmNlbnRhZ2UuCgpCRUQgZm9ybWF0IGZpbGUgbXVzdCBoYXZlIGEgaGVhZGVyIHJvdywgYW5kIHRoZSBjb29yZGluYXRlIHJlZmVyZW5jZSBmb3IgdGhlIGxvY3VzIGlzIDAuCgpgcmVhZEJlZEZpbGVgIGlzIGEgZnVuY3Rpb24gdGhhdCByZXNpZGVzIGluIHRoZSBuYW1lc3BhY2Ugb2YgTTNELiBUaGUgZm9sbG93aW5nIGlzIGl0cyBjb2RlLgoKYGBgCmZ1bmN0aW9uIChmaWxlcywgY29sRGF0YSwgYmVkX3R5cGUgPSAiRW5jb2RlIiwgZURhdGEgPSBOYU4pIAp7CiAgICBpZiAobnJvdyhjb2xEYXRhKSAhPSBsZW5ndGgoZmlsZXMpKSB7CiAgICAgICAgc3RvcCgiUm93IG51bWJlciBvZiBjb2xEYXRhIG11c3QgZXF1YWwgbGVuZ3RoIG9mIGZpbGVzLiIpCiAgICB9CiAgICBpZiAoIXRlc3RfY29sRGF0YV9zdHJ1Y3R1cmUoY29sRGF0YSkpIHsKICAgICAgICBzdG9wKCJGaWxlcyBtdXN0IGJlIHNwZWNpZmllZCBncm91cCBieSBncm91cC4gUGxlYXNlIHJlLW1ha2UgY29sRGF0YSB0b1xuICAgICAgICAgICAgIHNob3cgYWxsIG9mIG9uZSBncm91cCwgdGhlbiBhbGwgb2YgdGhlIG5leHQsIGV0Yy4iKQogICAgfQogICAgbWVzc2FnZSgiUmVhZGluZyBpbiBmaWxlc1xuIikKICAgIG1ldGhEYXRhIDwtIHJlYWRCZWRSYXcoZmlsZXMsIGNvbERhdGEsIGJlZF90eXBlID0gYmVkX3R5cGUpCiAgICBjYXQoIkJ1aWxkaW5nIEJTcmF3IG9iamVjdC5cbiIpCiAgICBmRGF0YSA8LSBtZXRoRGF0YVtbMV1dCiAgICBpZiAobGVuZ3RoKG1ldGhEYXRhKSA+IDEpIHsKICAgICAgICBmb3IgKGkgaW4gc2VxKGFsb25nID0gbWV0aERhdGEpWy0xXSkgewogICAgICAgICAgICBmRGF0YSA8LSB1bmlxdWUoYyhmRGF0YSwgbWV0aERhdGFbW2ldXSkpCiAgICAgICAgfQogICAgfQogICAgZWxlbWVudE1ldGFkYXRhKGZEYXRhKSA8LSBOVUxMCiAgICBuYW1lcyhmRGF0YSkgPC0gYXMuY2hhcmFjdGVyKDE6bGVuZ3RoKGZEYXRhKSkKICAgIHRSZWFkcyA8LSBtYXRyaXgoaW50ZWdlcihsZW5ndGggPSBsZW5ndGgoZkRhdGEpICogbGVuZ3RoKG1ldGhEYXRhKSksIAogICAgICAgIG5yb3cgPSBsZW5ndGgoZkRhdGEpKQogICAgbVJlYWRzIDwtIG1hdHJpeChpbnRlZ2VyKGxlbmd0aCA9IGxlbmd0aChmRGF0YSkgKiBsZW5ndGgobWV0aERhdGEpKSwgCiAgICAgICAgbnJvdyA9IGxlbmd0aChmRGF0YSkpCiAgICBmb3IgKGkgaW4gc2VxKGFsb25nID0gbWV0aERhdGEpKSB7CiAgICAgICAgbXRjaCA8LSBmaW5kT3ZlcmxhcHMoZkRhdGEsIG1ldGhEYXRhW1tpXV0pCiAgICAgICAgbXRjaC5tIDwtIGFzLm1hdHJpeChtdGNoKQogICAgICAgIGluZDEgPC0gbXRjaC5tWywgMV0KICAgICAgICBpbmQyIDwtIG10Y2gubVssIDJdCiAgICAgICAgdFJlYWRzW2luZDEsIGldIDwtIGVsZW1lbnRNZXRhZGF0YShtZXRoRGF0YVtbaV1dKSRyZWFkc1tpbmQyXQogICAgICAgIG1SZWFkc1tpbmQxLCBpXSA8LSBlbGVtZW50TWV0YWRhdGEobWV0aERhdGFbW2ldXSkkbWV0aHlsYXRlZFtpbmQyXQogICAgfQogICAgY29sbmFtZXModFJlYWRzKSA8LSByb3duYW1lcyhjb2xEYXRhKQogICAgY29sbmFtZXMobVJlYWRzKSA8LSByb3duYW1lcyhjb2xEYXRhKQogICAgcm93bmFtZXModFJlYWRzKSA8LSBuYW1lcyhmRGF0YSkKICAgIHJvd25hbWVzKG1SZWFkcykgPC0gbmFtZXMoZkRhdGEpCiAgICBpZiAoaXMubmEoZURhdGEpKSB7CiAgICAgICAgZURhdGEgPC0gU2ltcGxlTGlzdChleHAgPSAiRGVtbyIpCiAgICB9CiAgICBycmJzID0gQlNyYXcoZXhwdERhdGEgPSBlRGF0YSwgY29sRGF0YSA9IGNvbERhdGEsIHJvd1JhbmdlcyA9IGZEYXRhLCAKICAgICAgICB0b3RhbFJlYWRzID0gdFJlYWRzLCBtZXRoUmVhZHMgPSBtUmVhZHMpCiAgICByZXR1cm4ocnJicykKfQo8Ynl0ZWNvZGU6IDB4NTVkZjMxN2JlNDEwPgo8ZW52aXJvbm1lbnQ6IG5hbWVzcGFjZTpNM0Q+CmBgYAoKTmV4dCwgdXNlIEJpU2VxIHBhY2thZ2UgdG8gZ2VuZXJhdGUgR1JhbmdlcyBmaWxlIGZvciB0aGUgcmVnaW9uLiBCaVNlcSBwYWNrYWdlIGlzIGJ1aWx0IGZvciBwcm9jZXNzaW5nIGJpc3VsZml0ZS1jb252ZXJzaW9uIHNlcXVlbmNlZCBkYXRhLgoKYGBge3J9CnJyYnMuQ3BHcyA8LSBjbHVzdGVyU2l0ZXMob2JqZWN0ID0gcnJicywgZ3JvdXBzID0gdW5pcXVlKGNvbERhdGEocnJicykkZ3JvdXApLCBwZXJjLnNhbXBsZXMgPSAzLzQsIG1pbi5zaXRlcyA9IDIwLCBtYXguZGlzdCA9IDEwMCkKYGBgCgpgYGB7cn0KQ3BHcyA8LSBjbHVzdGVyU2l0ZXNUb0dSKHJyYnMuQ3BHcykKYGBgCgpJbiBvcmRlciB0byB1bmRlcnN0YW5kIHRoZSBuZXh0IHN0ZXAsIHlvdSBnb3R0YSBrbm93IHdoYXQgKipjb3ZlcmFnZSoqIChvciAqKnJlYWQgZGVwdGgqKikgaXMuIENvdmVyYWdlIGluIEROQSBzZXF1ZW5jaW5nIGlzIHRoZSBudW1iZXIgb2YgdW5pcXVlIHJlYWRzIHRoYXQgaW5jbHVkZSBhIGdpdmVuIG51Y2xlb3RpZGUgaW4gdGhlIHJlY29uc3RydWN0ZWQgIHNlcXVlbmNlLiBEZWVwIHNlcXVlbmNpbmcgYWltcyBhdCBoaWdoZXIgY292ZXJhZ2Ugb2YgZWFjaCB1bmlxdWUgcmVhZHMuCgpOb3csIHdlIGN1dCBvZmYgdGhlIHJlZ2lvbnMgd2l0aCBjb3ZlcmFnZSAocmVhZCBkZXB0aCkgbGVzcyB0aGFuIDEwMC4gYHRvdGFsUmVhZHMocnJicy5DcEdzKWAgaXMgYSBtYXRyaXggdGhhdCBjb250YWlucyBwb3NpdGlvbiBhbmQgdG90YWwgY292ZXJhZ2Ugb2YgZWFjaCBwb3NpdGlvbi4KCkluIHRoZSBmb2xsb3dpbmcgY29kZSwgYG92ZXJsYXBzYCBpcwpgYGAKPiBvdmVybGFwcwpIaXRzIG9iamVjdCB3aXRoIDU3NjEwNiBoaXRzIGFuZCAwIG1ldGFkYXRhIGNvbHVtbnM6CiAgICAgICAgICAgcXVlcnlIaXRzIHN1YmplY3RIaXRzCiAgICAgICAgICAgPGludGVnZXI+ICAgPGludGVnZXI+CiAgICAgICBbMV0gICAgICAgICAxICAgICAgICAgICAxCiAgICAgICBbMl0gICAgICAgICAxICAgICAgICAgICAyCiAgICAgICBbM10gICAgICAgICAxICAgICAgICAgICAzCiAgICAgICBbNF0gICAgICAgICAxICAgICAgICAgICA0CiAgICAgICBbNV0gICAgICAgICAxICAgICAgICAgICA1CiAgICAgICAuLi4gICAgICAgLi4uICAgICAgICAgLi4uCiAgWzU3NjEwMl0gICAgIDE0NDUzICAgICAgNTc2MTAyCiAgWzU3NjEwM10gICAgIDE0NDUzICAgICAgNTc2MTAzCiAgWzU3NjEwNF0gICAgIDE0NDUzICAgICAgNTc2MTA0CiAgWzU3NjEwNV0gICAgIDE0NDUzICAgICAgNTc2MTA1CiAgWzU3NjEwNl0gICAgIDE0NDUzICAgICAgNTc2MTA2CiAgLS0tLS0tLQpgYGAKYW5kIApgYGAKc3ViamVjdEhpdHMob3ZlcmxhcHNbIHF1ZXJ5SGl0cyhvdmVybGFwcyk9PSAxXSkKIFsxXSAgMSAgMiAgMyAgNCAgNSAgNiAgNyAgOCAgOSAxMCAxMSAxMiAxMyAxNCAxNSAxNiAxNyAxOCAxOSAyMCAyMSAyMiAyMyAyNCAyNSAyNiAyNyAyOCAyOSAzMCAzMSAzMiAzMyAzNCAzNSAzNgpgYGAKYHN1YmplY3RIaXRzYCBicmluZ3MgaW5kaWNlcyBvZiB3aGVyZSB0aGF0IG51bWJlciBvZiBvdmVybGFwIGlzLgoKYGBge3J9CmluZHMgPC0gdmVjdG9yKCkKb3ZlcmxhcHMgPC0gZmluZE92ZXJsYXBzKENwR3MsIHJvd1JhbmdlcyhycmJzLkNwR3MpKQpmb3IgKGkgaW4gMTpsZW5ndGgoQ3BHcykpewogIGNvdnMgPC0gY29sU3Vtcyh0b3RhbFJlYWRzKHJyYnMuQ3BHcylbc3ViamVjdEhpdHMob3ZlcmxhcHNbIHF1ZXJ5SGl0cyhvdmVybGFwcykgPT0gaV0pLF0pCiAgaWYgKCFhbnkoY292cyA8PSAxMDApKXsKICAgIGluZHMgPC0gYyhpbmRzLCBpKQogIH0KfQpDcEdzIDwtIENwR3NbaW5kc10Kcm0oaW5kcykKYGBgCgpGb3IgbmV4dCBhbmFseXNpcywgdXNlIG9ubHkgMTAwMCByZWdpb25zLiBVcGRhdGUgYG92ZXJsYXBzYCBvYmplY3QgYW5kIGV4cGxhaW4gY29ycmVjdCBgb3ZlcmxhcHNgLgoKYGBge3J9CnJyYnMgPC0gcnJicy5DcEdzCkNwR3MgPC0gQ3BHc1sxOjEwMDBdCm92ZXJsYXBzIDwtIGZpbmRPdmVybGFwcyhDcEdzLCByb3dSYW5nZXMocnJicykpCnJyYnMgPC0gcnJic1tzdWJqZWN0SGl0cyhvdmVybGFwcyldCm92ZXJsYXBzIDwtIGZpbmRPdmVybGFwcyhDcEdzLCByb3dSYW5nZXMocnJicykpCmBgYAoKClVzZSB0aGUgdHdvIGRhdGEgaW4gRU5DT0RFIGNvbnNvcnRpdW0gYEgxLWhFU0NgIGFuZCBgSzU2MmAgYXMgdGhlIHR3byBSUkJTIGZpbGVzLgoKVXNlIEJpU2VxIHBhY2thZ2UncyBmdW5jdGlvbnMgYGNsdXN0ZXJTaXRlc2AgYW5kIGBjbHVzdGVyU2l0ZXNUb0dSYCBmb3IgY2x1c3RlcmluZy4KCkRlbGV0ZSBpc2xhbmRzIG9mIGFsbCBzYW1wbGVzIHRoYXQgaGF2ZSBsZXNzIHRoYW4gMTAwPyBhbmQgaW5jbHVkZSBvbmx5IHRoZSBmaXJzdCB0aG91c2FuZC4KCkFuYWx5c2lzIGlzIHN0b3JlZCBpbiBHUmFuZ2VzIGFuZCBlYWNoIGl0ZW0gc2hvd3MgYSByZWdpb24ncyBzdGFydCBhbmQgZW5kLgoKYGBge3J9CnJyYnMKYGBgCgpgYGB7cn0KQ3BHcwpgYGAKIyMgTWF4aW11bSBNZWFuIERpc2NyZXBhbmN5IChNTUQpIGFuZCBDb3ZlcmFnZQoKTTNEIGFuYWx5c2lzIHVzZXMgTU1EIHZhbHVlcyBmb3IgdGhlIGVudGlyZSBtZXRoeWxhdGlvbiBkYXRhIGFuZCB0aGUgYXBwbGllZCByZWdpb25zLiBCb3RoIGNhbiBiZSBvYnRhaW5lZCB0aHJvdWdoIGBNM0RfV3JhcHBlcmAgZnVuY3Rpb24uIFNhbXBsZSBwYWlycycgZW50aXJlIG1ldGh5bGF0ZWQgTU1EIG1hdHJpeCBpcyB0aGUgcmV0dXJuIHZhbHVlIG9mIGBNM0RfV3JhcHBlcmAuIFRoZSBzZWNvbmQgcmVzdWx0IGlzIGVxdWl2YWxlbnQgY292ZXJhZ2UuIGBNTURsaXN0JEZ1bGxgIGFuZCBgTU1EbGlzdCRDb3ZlcmFnZWAgdmFyaWFibGVzIGFyZSB0aG9zZS4gU3Vic3RpdHV0ZSB0aGUgZGlmZmVyZW5jZSBiZXR3ZWVuIHRoZSB0d28gdmFsdWVzIGludG8gTTNEc3RhdCwgYW5kIHRoaXMgc2hvd3Mgc3RydWN0dXJlIG9mIGVhY2ggaXRlbXMuCgpgYGB7ciwgaW5jbHVkZT1GQUxTRX0KTU1EbGlzdCA8LSBNM0RfV3JhcHBlcihycmJzLCBvdmVybGFwcykKYGBgCgpgYGB7cn0KaGVhZChNTURsaXN0JEZ1bGwpCmBgYAoKYGBge3J9CmhlYWQoTU1EbGlzdCRDb3ZlcmFnZSkKYGBgCgpgYGB7cn0KTTNEc3RhdCA8LSBNTURsaXN0JEZ1bGwgLSBNTURsaXN0JENvdmVyYWdlCmBgYAoKQmVsb3cgaXMgYSBzY2F0dGVycGxvdCBmb3IgY29tcGFyaXNvbiBhbW9uZyByZXBsaWNhdGVzIG9mIEZ1bGwgTU1EIGFuZCBDb3ZlcmFnZSBNTUQuCmBgYHtyfQpyZXBDb2xzIDwtIGMoMSwgNikKcGxvdChhcy52ZWN0b3IoTU1EbGlzdCRGdWxsWywgcmVwQ29sc10pLCBhcy52ZWN0b3IoTU1EbGlzdCRDb3ZlcmFnZVssIHJlcENvbHNdKSwgeGxhYiA9ICdGdWxsIE1NRCcsIHlsYWIgPSAnQ292ZXJhZ2UgTU1EJywgbWFpbiA9ICdUZXN0IFN0YXRpc3RpYzogUmVwbGljYXRlIENvbXBhcmlzb24nKQphYmxpbmUoYSA9IDAsIGIgPSAxLCBjb2w9J3JlZCcsIGx3ZD0yKQpgYGAKYGBge3J9Cmdyb3VwQ29scyA8LSAyOjUKcGxvdChhcy52ZWN0b3IoTU1EbGlzdCRGdWxsWywgZ3JvdXBDb2xzXSksIGFzLnZlY3RvcihNTURsaXN0JENvdmVyYWdlWywgZ3JvdXBDb2xzXSksIHhsYWIgPSAnRnVsbCBNTUQnLCB5bGFiID0gJ0NvdmVyYWdlIE1NRCcsIG1haW4gPSAnVGVzdCBTdGF0aXN0aWNzOiBJbnRlci1Hcm91cCBDb21wYXJpc29uJykKYWJsaW5lKGE9MCwgYj0xLCBjb2w9J3JlZCcsIGx3ZD0yKQpgYGAKVGhlIGZvbGxvd2luZyBjb2RlIHByaW50cyBNM0QgYW5hbHlzaXMgcmVzdWx0IGFzIGhpc3RvZ3JhbXMuCmBgYHtyfQpyZXBDb2xzIDwtIGMoMTo2KQpncm91cENvbHMgPC0gMjo1Ck0zRHN0YXQgPC0gTU1EbGlzdCRGdWxsIC0gTU1EbGlzdCRDb3ZlcmFnZQpicmVha3MgPC0gc2VxKC0wLjIsIDEuMiwgMC4xKQpoUmVwcyA8LSBoaXN0KE0zRHN0YXRbLCByZXBDb2xzXSwgYnJlYWtzID0gYnJlYWtzLCBwbG90ID0gRkFMU0UpCmhHcm91cHMgPC0gaGlzdChyb3dNZWFucyhNM0RzdGF0WywgZ3JvdXBDb2xzXSksIGJyZWFrcyA9IGJyZWFrcywgcGxvdCA9IEZBTFNFKQpjb2wxIDwtICIjMDAwMEZGNDAiCmNvbDIgPC0gIiNGRjAwMDA0MCIKcGxvdChoUmVwcywgY29sID0gY29sMSwgZnJlcSA9IEZBTFNFLCB5bGFiID0gIkRlbnNpdHkiLCB4bGFiID0gIk1NRCBzdGF0aXN0aWMiLCBtYWluID0gIk0zRCBTdGF0IEhpc3RvZ3JhbSIpCnBsb3QoaEdyb3VwcywgY29sID0gY29sMiwgZnJlcSA9IEZBTFNFLCBhZGQgPSBUUlVFKQpsZWdlbmQoeCA9IDAuNSwgeSA9IDMsIGxlZ2VuZCA9IGMoIlJlcGxpY2F0ZXMiLCAiR3JvdXBzIiksIGZpbGwgPSBjKGNvbDEsIGNvbDIpKQpgYGAKCk5leHQgaXMgaGlzdG9uZSBtb2RpZmljYXRpb24gZGF0YSBzZWFyY2ggaW4gdGhlIEVwaUZhY3RvcnMgZGF0YWJhc2Uu